Código
install.packages(c(
"sf", "spdata", "RColorBrewer", "spdep", "spatialreg",
"readxl", "dplyr", "ggplot2", "geobr", "ggspatial", "ggthemes",
"ggrepel", "viridis", "tmap", "car", "janitor"
))Dados de áreas (areal data) são observações agregadas em unidades espaciais discretas, como bairros, municípios ou regiões administrativas. Essas unidades são representadas por polígonos em um sistema de informação geográfica (SIG). Exemplos comuns incluem taxas de criminalidade por bairro, incidência de doenças por município e dados socioeconômicos por região.
As áreas podem ser regulares ou irregulares. Em geral, áreas representadas por regiões geográficas como estados e municípios são denominadas de espacialmente irregulares, enquanto áreas representadas por pixels, como imagens de sensoriamento remoto, são denominadas de espacialmente regulares (formado por uma grade de espaços iguais) (Scalon 2024).
Nesta abordagem, modelamos uma variável resposta agregada em polígonos (ex: taxas de crime por bairro) usando variáveis independentes, corrigindo a autocorrelação espacial via matriz de pesos \(W\).
Os dados das variáveis aleatórias podem ser de qualquer tipo (contínuas, discretas, etc) e tipicamente, representam toda a área, ou seja, não se dispõe da localização exata dos eventos dentro da área, mas sim de um valor que foi agregado para aquela área. Embora esses dados estejam associados com toda a área, em geral, eles são atribuídos a um ponto específico dentro da área, esse ponto pode ser o centroide da área (Scalon 2024).
O objetivo da análise de dados de áreas, é o mesmo da análise dos outros tipos de dados espaciais, ou seja, caracterizar o processo estocástico que gerou os dados (Scalon 2024).
Na análise de dados de áreas, como existe apenas um valor \(z_i\) para a área \(W_i\), deve-se assumir algum tipo de estacionariedade e para alguns casos será necessário assumir a normalidade multivariada (Scalon 2024).
Para esse tipo de análise de dados espaciais, dados de áreas, pode estar interessado em um ou mais objetivos. Por exemplo, Scalon (2024) cita alguns dos objetivos, que podem ser:
- Explorar o comportamento das variáveis aleatórias na região de estudo; - Obter mapas mais suaves do que o mapa dos dados observados; - etc…
Para a realização de estatística espacial, é necessário que os dados em análise possuam o termo de autocorrelação espacial entre a variável aleatória.
O software para este minicurso, será o R, usando a IDE Rstudio.
Formato vetorial poligonal, como shapefiles (.shp) ou arquivos GeoPackage (.gpkg), que armazenam as geometrias dos polígonos e seus atributos associados.
Para realizar a análise de dados de áreas, utilizaremos os seguintes pacotes:
Execute o código abaixo para instalar as bibliotecas:
sf: Para manipulação de dados espaciais (Simple Features).spdata: Contém conjuntos de dados espaciais para exemplos.RColorBrewer: Para criação de paletas de cores em mapas.spdep: Para criar matrizes de pesos e testes de autocorrelação.spatialreg: Para ajustar modelos espaciais autoregressivos.readxl: Para leitura de arquivos Excel (.xls, .xlsx).dplyr: Para manipulação de dados.ggplot2: Para criação de gráficos e visualização de dados.geobr: Para baixar shapefiles oficiais do Brasil (IBGE).ggspatial: Adiciona elementos espaciais (norte, escala) ao ggplot2.ggthemes: Temas adicionais para gráficos.ggrepel: Para rótulos de texto que não se sobrepõem.viridis: Paletas de cores acessíveis para visualização de dados.tmap: Para criação de mapas temáticos flexíveis.car: Para diagnósticos de regressão (ex: VIF, QQ-plots).janitor: Para limpeza de dados (ex: nomes de colunas).Obtendo ajuda: Para consultar a documentação de um pacote ou função, utilize os comandos help() ou ?.
O modelo espacial autoregressivo (SAR) é descrito conforme (Anselin1988?)
\[\begin{equation} Y = \delta WY + X\beta + \epsilon \label{eq:sar} \end{equation}\]
Onde
\(Y\) é um vetor de observações ( n \(\times\) 1 ) nas n áreas
\(X\) é a matriz das variáveis independentes, de tamanho ( n \(\times\) i ),
\(\delta\) é o coeficiente de defasagem da variável dependente espacial,
\(\beta\) é o vetor de parâmetros da regressão, de tamanho ( i \(\times\) 1 ),
\(W\) é a matriz de ponderação, de tamanho ( n \(\times\) n ),
\(\epsilon\) é o vetor de erros, de tamanho ( n \(\times\) 1 ) não correlacionados que seguem uma distribuição normal com média zero e variância constante, isto é, \(\epsilon \sim N(0, I\sigma^2\))
Abaixo apresentamos o passo a passo da análise, utilizando dados reais de incidência.
Primeiramente, carregamos todas as bibliotecas que serão utilizadas na análise.
# Carregando dados de incidência
dados <- read_excel("C:/Meus_documentos/Congresso_MGEST_2025/Chikugunha/dados/chiku.xlsx")
dados <- janitor::clean_names(dados)
dados <- dados |>
select(cod, muni, x2023_0) |>
mutate(x2023_0 = as.numeric(replace(x2023_0, x2023_0 == "-", 0)))
# Carregando dados de densidade/população
dados1 <- read_excel("C:/Meus_documentos/Congresso_MGEST_2025/Chikugunha/dados/densidad.xlsx", sheet = 3)
dados1 <- dados1 |>
separate(muni, into = c("Municipio", "Estado"), sep = " \\(") %>%
mutate(Estado = gsub("\\)", "", Estado))
dados1 <- dados1[-c(1, 219:222),-2 ] # Limpeza de linhas# Baixar malha municipal (MA)
todos <- read_municipality(code_muni = "all", year = 2022, showProgress = FALSE)
ma <- filter(todos, abbrev_state == "MA")
ma <- ma |>
mutate(code_muni = str_sub(code_muni, 1, -2),
code_muni = as.numeric(code_muni))
# Unificação (Join)
dados_juntos <- ma |>
left_join(dados, by = c("code_muni" = "cod"))
dados_juntos <- dados_juntos[, -c(4,5,7,8)]
dad <- cbind(dados_juntos, dados1) |>
mutate(across(everything(), ~replace_na(., 0))) |>
mutate(popu = as.numeric(popu))
dad <- dad[, -c(3,4,6)]
# Cálculo da taxa de incidência
dad <- dad |>
mutate(inc_100 = (x2023_0 / popu)*10000, .after =x2023_0,
log_inc_100 = log((x2023_0 / popu)*10000))# Gráficos exploratórios
par(mfrow = c(1,2))
boxplot(dad$inc_100, main = "Incidência")
hist(dad$inc_100, main = "Histograma")
# Regressão Linear (OLS)
modelo_lm_2 <- lm(inc_100 ~ area_da_unidade + densidade_demográfica, data = dad)
summary(modelo_lm_2)
# Diagnóstico
vif(modelo_lm_2)
qqPlot(modelo_lm_2$residuals, main = "QQ Plot OLS")# Matriz de vizinhança e pesos
nb <- poly2nb(dad$geom, queen=TRUE)
lw <- nb2listw(nb, style="W", zero.policy=TRUE)
# Teste de Moran Global
moran.test(dad$inc_100, lw, randomisation = TRUE)
# Diagrama de Dispersão de Moran
moran.plot(modelo_lm_2$residual, lw, main = "Diagrama de Dispersão de Moran (OLS)")# Modelo SAR (Spatial Lag)
modelo_sar1 <- lagsarlm(inc_100 ~ area_da_unidade + densidade_demográfica,
data = dad, listw = lw, zero.policy = TRUE)
summary(modelo_sar1)
# Modelo SEM (Spatial Error) - Referido como CAR no código original
modelo_car1 <- errorsarlm(inc_100 ~ area_da_unidade + densidade_demográfica,
data = dad, listw = lw, zero.policy = TRUE)
summary(modelo_car1)# Adicionando resultados ao mapa
dad$pred_sar <- modelo_sar1$fitted.values
dad$resid_sar <- residuals(modelo_sar1)
# Mapa de Predição
ggplot(dad) +
geom_sf(aes(fill = pred_sar), color = NA) +
scale_fill_viridis_c(option = "magma", name = "Predição SAR") +
annotation_scale(location = "bl", width_hint = 0.4) +
annotation_north_arrow(location = "tl", which_north = "true") +
theme_minimal()
# Clusters LISA
dad_sf <- st_as_sf(dad)
lisa <- localmoran(dad_sf$inc_100, lw)
dad_sf$lisa_i <- lisa[,1]
dad_sf$pval <- lisa[,5]
dad_sf <- dad_sf %>%
mutate(cluster = case_when(
lisa_i > 0 & pval <= 0.05 & inc_100 > mean(inc_100) ~ "Alto-Alto",
lisa_i > 0 & pval <= 0.05 & inc_100 < mean(inc_100) ~ "Baixo-Baixo",
lisa_i < 0 & pval <= 0.05 & inc_100 > mean(inc_100) ~ "Alto-Baixo",
lisa_i < 0 & pval <= 0.05 & inc_100 < mean(inc_100) ~ "Baixo-Alto",
TRUE ~ "Não significativo"
))
tm_shape(dad_sf) +
tm_polygons("cluster", palette = c("red", "blue", "white", "lightblue", "gray80"))